<!DOCTYPE html PUBLIC "-//W3C//DTD XHTML 1.0 Transitional//EN"
  "http://www.w3.org/TR/xhtml1/DTD/xhtml1-transitional.dtd">


<html xmlns="http://www.w3.org/1999/xhtml">
  <head>
    <meta http-equiv="Content-Type" content="text/html; charset=utf-8" />
    
    <title>第14章 使用Bio.motifs进行模体序列分析 &mdash; Biopython_en 1.0 documentation</title>
    
    <link rel="stylesheet" href="_static/default.css" type="text/css" />
    <link rel="stylesheet" href="_static/pygments.css" type="text/css" />
    
    <script type="text/javascript">
      var DOCUMENTATION_OPTIONS = {
        URL_ROOT:    './',
        VERSION:     '1.0',
        COLLAPSE_INDEX: false,
        FILE_SUFFIX: '.html',
        HAS_SOURCE:  true
      };
    </script>
    <script type="text/javascript" src="_static/jquery.js"></script>
    <script type="text/javascript" src="_static/underscore.js"></script>
    <script type="text/javascript" src="_static/doctools.js"></script>
    <link rel="top" title="Biopython_en 1.0 documentation" href="index.html" />
    <link rel="next" title="Chapter 15 Cluster analysis" href="chr15.html" />
    <link rel="prev" title="Chapter 13 Phylogenetics with Bio.Phylo" href="chr13.html" /> 
  </head>
  <body>
    <div class="related">
      <h3>Navigation</h3>
      <ul>
        <li class="right" style="margin-right: 10px">
          <a href="genindex.html" title="General Index"
             accesskey="I">index</a></li>
        <li class="right" >
          <a href="chr15.html" title="Chapter 15 Cluster analysis"
             accesskey="N">next</a> |</li>
        <li class="right" >
          <a href="chr13.html" title="Chapter 13 Phylogenetics with Bio.Phylo"
             accesskey="P">previous</a> |</li>
        <li><a href="index.html">Biopython_en 1.0 documentation</a> &raquo;</li> 
      </ul>
    </div>  

    <div class="document">
      <div class="documentwrapper">
        <div class="bodywrapper">
          <div class="body">
            
  <div class="section" id="bio-motifs">
<h1>第14章   使用Bio.motifs进行模体序列分析<a class="headerlink" href="#bio-motifs" title="Permalink to this headline">¶</a></h1>
<p>这章为Biopython中的``Bio.motifs`` 包进行一个主要的介绍。这个包是为了方便那些需要进行模体序列分析的人们而特意提供的，所以我想你们在使用时肯定对模体序列分析的一些相关要点都很熟悉。假如在使用中遇到不清楚的地方，请您查阅 <a class="reference external" href="#sec:links">14.8</a> 相关章节以获得有关的信息。</p>
<p>这章的大部分内容是介绍Biopython 1.61 之前版本中新加入的 <tt class="docutils literal"><span class="pre">Bio.motifs</span></tt> 包，该包替代了Biopython 1.50版本中的 <tt class="docutils literal"><span class="pre">Bio.Motif</span></tt> 包，而``Bio.Motif`` 包是基于较早版本的Biopython 中的两个模块 <tt class="docutils literal"><span class="pre">Bio.AlignAce</span></tt> 和 <tt class="docutils literal"><span class="pre">Bio.MEME</span></tt> 。``Bio.motifs`` 包较好地综合了上述的几个模块的功能，做为一个统一模块工具。</p>
<p>说到其他库，看到这里，你或许会对 <a class="reference external" href="http://fraenkel.mit.edu/TAMO/">TAMO</a> 感兴趣，这是另一个分析模体序列的Python库。它能提供更多关于*de-novo*模体的查找方式，不过它并没有纳入到Biopython中，而且在商业用途上还有一些限制。</p>
<div class="section" id="id1">
<h2>14.1  模体对象<a class="headerlink" href="#id1" title="Permalink to this headline">¶</a></h2>
<p>由于我们感兴趣的是模体分析，所以我们需要先看看 <tt class="docutils literal"><span class="pre">Motif</span></tt> 对象。对此我们需要先导入Bio.motifs库：</p>
<p>然后我们可以开始创建我们第一个模体对象。我们可以从模体的实例列表中创建一个 <tt class="docutils literal"><span class="pre">Motif</span></tt> 对象，也可以通过读取模体数据库中或模体查找软件产生的文件来获得一个 <tt class="docutils literal"><span class="pre">Motif</span></tt> 对象。</p>
<div class="section" id="id2">
<h3>14.1.1  从实例中创建一个模体<a class="headerlink" href="#id2" title="Permalink to this headline">¶</a></h3>
<p>假设我们有一些DNA模体的实例：</p>
<p>然后我们可以如下创建一个模体对象：</p>
<p>这些实例被存储在一个名为 <tt class="docutils literal"><span class="pre">m.instances</span></tt> 的属性中，这个其实也就是一个Python的列表，只不过附加了一些功能，这些功能将在之后介绍。将这些模体对象打印出来后就可以看出这些实例是从哪构建出来的。</p>
<p>模体的长度像其他一些实例一些被定义为序列的长度：</p>
<p>模体对象有一个 <tt class="docutils literal"><span class="pre">.counts</span></tt> 属性，可以用来查看碱基在每个位置的数目。可以把这个统计表用易读的格式打印出来：</p>
<p>你也可以像使用字典一样获取这些数目：</p>
<p>但是你也可以把它看成一个二维数列，核苷酸作为列，位置作为行：</p>
<p>你还可以直接获得核苷酸数目矩阵中的列</p>
<p>除了使用核苷酸本身，你还可以使用模体碱基序列按字符排序后的核苷酸索引：</p>
<p>模体有一个相关联的一致序列，这个序列被定义为由 <tt class="docutils literal"><span class="pre">.counts</span></tt> 矩阵相应列中具有最大值的碱基，这些碱基是按模体序列排列的：</p>
<p>反一致序列也一样，只不过是由 <tt class="docutils literal"><span class="pre">.counts</span></tt> 矩阵中相应列的最小值来选：</p>
<p>你也可以利用降级的一致序列，用不确定核苷酸来表示序列某一位置的所有核苷酸：</p>
<p>在这，W和R都是按照IUPAC不确定核苷酸表规定的：W代表A或T，V代表A，C或G [<a class="reference external" href="#cornish1985">10</a>] 。这些降级一致序列是按照Cavener指定的规则[<a class="reference external" href="#cavener1987">11</a>]来建立的。</p>
<p>反向互补序列和降级一致序列都只在DNA模体中有。</p>
</div>
<div class="section" id="id3">
<h3>14.1.2  读取模体<a class="headerlink" href="#id3" title="Permalink to this headline">¶</a></h3>
<p>从实例手动创建一个模体确实有点无趣，所以用一些I/O函数来读写模体是很有用的。目前对于如何存储模体还没有一些真正的标准，不过有一些格式用得比其他一些经常。这其中最重要的区别在于模体表示是基于实例还是某种PWM矩阵。</p>
<div class="section" id="jaspar">
<h4>JASPAR<a class="headerlink" href="#jaspar" title="Permalink to this headline">¶</a></h4>
<p>作为一个最流行的模体数据库 <a class="reference external" href="http://jaspar.genereg.net">JASPAR</a> 它不是以一系列的实例就是频率矩阵。比如，下面就是JASPAR <tt class="docutils literal"><span class="pre">Arnt.sites</span></tt> 文件的开头和结尾行显示了老鼠螺旋-环-螺旋转录因子Arnt的结合位点：</p>
<p>那些用大字字母表示的序列的一部分就是被开来互相比对的模体实例。</p>
<p>我们可以从下面的实例创建一个 <tt class="docutils literal"><span class="pre">Motif</span></tt> 对象：</p>
<p>从这个模体创建的实例存储在该模体的 <tt class="docutils literal"><span class="pre">.instances</span></tt> 属性：</p>
<p>这个模体的计数矩阵可以自动的由这些实例中计算出来：</p>
<p>JASPAR数据库也可以让模体像计数矩阵一样获得，不需要那些创建它们的实例。比如，下面这个JASPAR文件 <tt class="docutils literal"><span class="pre">SRF.pfm</span></tt> 包含了人类SRF转录因子的计数矩阵：</p>
<p>我们可以像下面那样为计数矩阵创建一个模体：</p>
<p>由于这个模体是由计数矩阵直接创建的，所以它没有相关和实例：</p>
<p>我们可以获得这两个模体的一致序列：</p>
</div>
<div class="section" id="meme">
<h4>MEME<a class="headerlink" href="#meme" title="Permalink to this headline">¶</a></h4>
<p>MEME [<a class="reference external" href="#bailey1994">12</a>] 是一个用来在一堆相关DNA或蛋白质序列中发现模体的工具。它输入一组相关DNA或蛋白质序列，输出所要求的模体。因此和JASPAR文件相比，MEME输出文件里面一般是含有多个模体。以下是一个例子。</p>
<p>在输出文件的开头，有一些MEME生成的关于MEME和所用MEME版本的背景信息：</p>
<p>往下，简要概括了输入的训练序列集：</p>
<p>以及所使用到的命令：</p>
<p>接下来就是每个被发现模体的详细信息：</p>
<p>读取这个文件（以 <tt class="docutils literal"><span class="pre">meme.dna.oops.txt</span></tt> 存储），使用下面的方法：</p>
<p><tt class="docutils literal"><span class="pre">motifs.parse</span></tt> 命令直接读取整个文件，所以在使用后可以关闭这个文件。其中头文件信息被存储于属性中</p>
<p>这个数据记录是 <tt class="docutils literal"><span class="pre">Bio.motifs.meme.Record</span></tt> 类的一个对象。这个类继承于列表（list），所以你可以把这个 <tt class="docutils literal"><span class="pre">record</span></tt> 看成模体对象的一个列表：</p>
<p>除了一般的模体属性外，每个模体同时还保存着它们由MEME计算的各自特异信息。例如：</p>
<p>除了像上面所做的用索引来获得相关记录，你也可以用它的名称来找到这个记录：</p>
<p>每个模体都有一个 <tt class="docutils literal"><span class="pre">.instances</span></tt> 拥有创建这个模体序列实例的属性，能够为每个实例提供一些信息：</p>
</div>
<div class="section" id="mast">
<h4>MAST<a class="headerlink" href="#mast" title="Permalink to this headline">¶</a></h4>
</div>
<div class="section" id="transfac">
<h4>TRANSFAC<a class="headerlink" href="#transfac" title="Permalink to this headline">¶</a></h4>
<p>TRANSFAC是一个为转录因子手动创建的一个专业数据库，同时还包括染色体结合位点和DNA结合的描述 [<a class="reference external" href="#matys2003">27</a>]。TRANSFAC数据库中所用的文件格式在今天还被其他地方所用到，我们下面将介绍TRANSFAC文件格式。</p>
<p>TRANSFAC文件格式简单概括如下：</p>
<p>这个文件显示了模体 <tt class="docutils literal"><span class="pre">motif1</span></tt> 中12个核苷酸的频率矩阵。总的来说，一个TRANSFAC文件里面可以包含多个模体。以下是示例文件 <tt class="docutils literal"><span class="pre">transfac.dat</span></tt> 的内容：</p>
<p>可用下面的方法读取TRANSFAC文件：</p>
<p>如果有总版本号的话，它是存储在 <tt class="docutils literal"><span class="pre">record.version</span></tt> 中：</p>
<p>每个在 <tt class="docutils literal"><span class="pre">record</span></tt> 中的模体都是 <tt class="docutils literal"><span class="pre">Bio.motifs.transfac.Motif</span></tt> 类的实例，这些实例同时继承 <tt class="docutils literal"><span class="pre">Bio.motifs.Motif</span></tt> 类和Python字典的特性。这些字典用双字母的键来存储关于这个模体的其他附加信息：</p>
<p>TRANSFAC文件一般比这些例子更好的详细，包含了许多关于模体的附加信息。表格 <a class="reference external" href="#table:transfaccodes">14.1.2</a> 列出了在TRANSFAC文件常见的双字母含义：</p>
<blockquote>
<div><table border="1" class="docutils">
<colgroup>
<col width="100%" />
</colgroup>
<tbody valign="top">
<tr class="row-odd"><td>Table 14.1: TRANSFAC文件中常见的字段</td>
</tr>
</tbody>
</table>
<table border="1" class="docutils">
<colgroup>
<col width="16%" />
<col width="84%" />
</colgroup>
<tbody valign="top">
<tr class="row-odd"><td><tt class="docutils literal"><span class="pre">AC</span></tt></td>
<td>Accession numbers 序列号</td>
</tr>
<tr class="row-even"><td><tt class="docutils literal"><span class="pre">AS</span></tt></td>
<td>Accession numbers, secondary 第二序列号</td>
</tr>
<tr class="row-odd"><td><tt class="docutils literal"><span class="pre">BA</span></tt></td>
<td>Statistical basis 统计依据</td>
</tr>
<tr class="row-even"><td><tt class="docutils literal"><span class="pre">BF</span></tt></td>
<td>Binding factors 结合因子</td>
</tr>
<tr class="row-odd"><td><tt class="docutils literal"><span class="pre">BS</span></tt></td>
<td>Factor binding sites underlying the matrix
基于矩阵的转录结合位点</td>
</tr>
<tr class="row-even"><td><tt class="docutils literal"><span class="pre">CC</span></tt></td>
<td>Comments 注解</td>
</tr>
<tr class="row-odd"><td><tt class="docutils literal"><span class="pre">CO</span></tt></td>
<td>Copyright notice 版权事项</td>
</tr>
<tr class="row-even"><td><tt class="docutils literal"><span class="pre">DE</span></tt></td>
<td>Short factor description 短因子说明</td>
</tr>
<tr class="row-odd"><td><tt class="docutils literal"><span class="pre">DR</span></tt></td>
<td>External databases 外部数据库</td>
</tr>
<tr class="row-even"><td><tt class="docutils literal"><span class="pre">DT</span></tt></td>
<td>Date created/updated 创建或更新日期</td>
</tr>
<tr class="row-odd"><td><tt class="docutils literal"><span class="pre">HC</span></tt></td>
<td>Subfamilies 亚家庭名称</td>
</tr>
<tr class="row-even"><td><tt class="docutils literal"><span class="pre">HP</span></tt></td>
<td>Superfamilies 超家庭名称</td>
</tr>
<tr class="row-odd"><td><tt class="docutils literal"><span class="pre">ID</span></tt></td>
<td>Identifier 身份证</td>
</tr>
<tr class="row-even"><td><tt class="docutils literal"><span class="pre">NA</span></tt></td>
<td>Name of the binding factor 结合因子的名称</td>
</tr>
<tr class="row-odd"><td><tt class="docutils literal"><span class="pre">OC</span></tt></td>
<td>Taxonomic classification 分类</td>
</tr>
<tr class="row-even"><td><tt class="docutils literal"><span class="pre">OS</span></tt></td>
<td>Species/Taxon 种类或分类</td>
</tr>
<tr class="row-odd"><td><tt class="docutils literal"><span class="pre">OV</span></tt></td>
<td>Older version 旧版本</td>
</tr>
<tr class="row-even"><td><tt class="docutils literal"><span class="pre">PV</span></tt></td>
<td>Preferred version 首先版本</td>
</tr>
<tr class="row-odd"><td><tt class="docutils literal"><span class="pre">TY</span></tt></td>
<td>Type 类型</td>
</tr>
<tr class="row-even"><td><tt class="docutils literal"><span class="pre">XX</span></tt></td>
<td>Empty line; these are not stored in the Record.
空白行;没在记录中存储的数据</td>
</tr>
</tbody>
</table>
</div></blockquote>
<p>Each motif also has an attribute <tt class="docutils literal"><span class="pre">.references</span></tt> containing the
references associated with the motif, using these two-letter keys:
每个模体同时也有一个包含与这个模体相关参考资料的 <tt class="docutils literal"><span class="pre">references</span></tt> 属性，用下面的双字母键来获得：</p>
<blockquote>
<div><table border="1" class="docutils">
<colgroup>
<col width="100%" />
</colgroup>
<tbody valign="top">
<tr class="row-odd"><td>Table 14.2: TRANSFAC文件中用来存储参考资料的字段</td>
</tr>
</tbody>
</table>
<table border="1" class="docutils">
<colgroup>
<col width="24%" />
<col width="76%" />
</colgroup>
<tbody valign="top">
<tr class="row-odd"><td><tt class="docutils literal"><span class="pre">RN</span></tt></td>
<td>Reference number 参考数目</td>
</tr>
<tr class="row-even"><td><tt class="docutils literal"><span class="pre">RA</span></tt></td>
<td>Reference authors 参考资料作者</td>
</tr>
<tr class="row-odd"><td><tt class="docutils literal"><span class="pre">RL</span></tt></td>
<td>Reference data 参考数据</td>
</tr>
<tr class="row-even"><td><tt class="docutils literal"><span class="pre">RT</span></tt></td>
<td>Reference title 参考标题</td>
</tr>
<tr class="row-odd"><td><tt class="docutils literal"><span class="pre">RX</span></tt></td>
<td>PubMed ID</td>
</tr>
</tbody>
</table>
</div></blockquote>
<p>将TRANSFAC文件按原来格式打印出来：</p>
<p>通过用字符串形式来截取输出并且保存在文件中，你可以按TRANSFAC的格式导出这些模体：</p>
</div>
</div>
<div class="section" id="id4">
<h3>14.1.3  模体写入<a class="headerlink" href="#id4" title="Permalink to this headline">¶</a></h3>
<p>说到导出，我们可以先看看导出函数。以JASPAR <tt class="docutils literal"><span class="pre">.pfm</span></tt> 格式导出模体文件，可以用：</p>
<p>用TRANSFAC类似的格式导出一个模体，可以用：</p>
<p>你可以用 <tt class="docutils literal"><span class="pre">motifs.write</span></tt> 来写入多个模体。这个函数在使用的时候不必担心这些模体来自于TRANSFAC文件。比如：</p>
</div>
<div class="section" id="id5">
<h3>14.1.4  绘制序列标识图<a class="headerlink" href="#id5" title="Permalink to this headline">¶</a></h3>
<p>如何能够联网，我们可以在该网址 <a class="reference external" href="http://weblogo.berkeley.edu">weblogo</a> 创建一个：</p>
<p>将得到的标识图存储成PNG格式。</p>
</div>
</div>
<div class="section" id="id6">
<h2>14.2  位置权重矩阵<a class="headerlink" href="#id6" title="Permalink to this headline">¶</a></h2>
<p>模体对象的 <tt class="docutils literal"><span class="pre">.counts</span></tt> 属性能够显示在序列上每个位置核苷酸出现的次数。我们可以把这矩阵除以序列中的实例数目来标准化这矩阵，得到每个核苷酸在序列位置上出现概率。我们把这概率看作位置权重矩阵。不过，要知道在字面上，这个术语也可以用来说明位置特异性得分矩阵，这个我们将会在下面讨论。</p>
<p>通常来说，伪计数（pseudocounts）在归一化之前都加到了每个位置中。这样可以避免在这序列上过度拟合位置权重矩阵以至趋向于模体的实例的有限数量，还可以避免概率为0。向每个位置的核苷酸添加一个固定的伪计数，可以为 <tt class="docutils literal"><span class="pre">pseudocounts</span></tt> 参数指定一个数值：</p>
<p>另外， <tt class="docutils literal"><span class="pre">pseudocounts</span></tt> 可以利用字典为每个核苷酸指定一个伪计数值。例如，由于在人类基因组中GC含量大概为40%,因此可以选择下面这些伪计数值：</p>
<p>位置权重矩阵有它自己的方法计算一致序列、反向一致序列和降级的一致序列（degenerate consensus sequences)：</p>
<p>应当注意到由于伪计数的原因，由位置仅重矩阵计算得到的降级一致序列和由模体中实例计算得到的降级一致序列有一点不同：</p>
<p>位置权重矩阵的反向互补矩阵可以直接用 <tt class="docutils literal"><span class="pre">pwm</span></tt> 计算出来：</p>
</div>
<div class="section" id="id7">
<h2>14.3  位置特异性得分矩阵<a class="headerlink" href="#id7" title="Permalink to this headline">¶</a></h2>
<p>使用背景分布和加入伪计数的PWM，很容易就能计算出log-odds比率，提供特定标记的log odds值，这值来自于在这个背景的模体。我们可以用在位置仅重矩阵中 <tt class="docutils literal"><span class="pre">.log-odds()</span></tt> 方法：</p>
<p>这时我们可以更经常看到特定标记和背景下的正值和负值。0.0意味着在模体和在背景中看到一个标记的可能性是一样的。</p>
<p>上面是假设A,C,G和T在背景中出现的概率是相同的。那在A,C,G和T出现概率不同的情况下，为了计算特定背景下的位置特异性得分矩阵，可以使用 <tt class="docutils literal"><span class="pre">background</span></tt> 参数。例如，在40%GC含量的背景下，可以用：</p>
<p>从PSSM中得到的最大和最小值被存储在 <tt class="docutils literal"><span class="pre">.max</span></tt> 和 <tt class="docutils literal"><span class="pre">.min</span></tt> 属性中：</p>
<p>在特定背景下计算平均值和标准方差可以使用 <tt class="docutils literal"><span class="pre">.mean</span></tt> 和 <tt class="docutils literal"><span class="pre">.std</span></tt> 方法。</p>
<p>如果没有指定特定的背景，就会使用一个统一的背景。因为同KL散度或相对熵的值相同，所以平均值就显得特别重要，并且它也是同背景相比的模体信息含量的测量方法。由于在Biopython中用以2为底的对数来计算log-odds值，信息含量的的单位是bit。</p>
<p><tt class="docutils literal"><span class="pre">.reverse_complement</span></tt>, <tt class="docutils literal"><span class="pre">.consensus</span></tt>, <tt class="docutils literal"><span class="pre">.anticonsensus</span></tt> 和``.degenerate_consensus`` 方法可以直接对PSSM使用。</p>
</div>
<div class="section" id="id8">
<h2>14.4  搜索实例<a class="headerlink" href="#id8" title="Permalink to this headline">¶</a></h2>
<p>模体最常用的功能就是在一些序列中的查找它的实例。在这节，我们会用以下一个人造序列：</p>
<div class="section" id="id9">
<h3>14.4.1  搜索准确匹配实例<a class="headerlink" href="#id9" title="Permalink to this headline">¶</a></h3>
<p>查找实例最简单的方法就是查找模体实例的准确匹配：</p>
<p>我们可获得反向互补序列（找到互补链的实例）：</p>
</div>
<div class="section" id="pssm">
<h3>14.4.2  用PSSM得分搜索匹配实例<a class="headerlink" href="#pssm" title="Permalink to this headline">¶</a></h3>
<p>在模体中很容易找出相应的位置,引起对模体的高log-odds值：</p>
<p>负值的位置是指在测试序列的反向链中找到的模体的实例，而且得力于Python的索引。在 <tt class="docutils literal"><span class="pre">pos</span></tt> 的模体实例可以用 <tt class="docutils literal"><span class="pre">test_seq[pos:pos+len(m)]</span></tt> 来定位，不管 <tt class="docutils literal"><span class="pre">pos</span></tt> 值是正还是负。</p>
<p>你可能注意到阀值参数，在这里随意地设为3.0。这里是 <em>log</em><sub>2</sub> ，所以我们现在开始寻找那些在模体中出现概率为背景中出现概率8倍序列。默认的阀值是0.0,在此阀值下，会把所有比背景中出现概率大的模体实例都找出来。</p>
<p>通常来说，上面是计算PSSM得分的最快方法。这些得分只能由前导链用 <tt class="docutils literal"><span class="pre">pssm.calculate</span></tt> 。为了得到互补链的PSSM值，你可以利用PSSM的互补矩阵：</p>
</div>
<div class="section" id="id10">
<h3>14.4.3  选择一个得分阀值<a class="headerlink" href="#id10" title="Permalink to this headline">¶</a></h3>
<p>如果不想那么随意设定一个阀值，你可以探究一下PSSM得分的分布。由于得分的空间分布随着模体长度而成倍增长，我们正用一个近似于一个给定精度值来使计算成本可管理：</p>
<p><tt class="docutils literal"><span class="pre">distribution</span></tt> 对象可以用来决定许多不同的阀值。我们可以指定一个所要求的假阳性率（找到一个由这个序列在这个背景下产生的一个模体实例的概率）：</p>
<p>或者假阴性率（找不到模体产生的实例概率）：</p>
<p>或者一个阀值（近似），满足假阳性率和假阴性率之间的关系（fnr/fpr≃ <em>t</em>)：</p>
<p>或者一个阀值能够大体满足假阳性率和信息含量的 −<em>log</em> 值之间的相等关系（像Hertz和Stormo的Patser软件所用的一样）：</p>
<p>比如在我们这个模体中，当以1000比率的平衡阀值查找实例，你可以得到一个让你获得相同结果的阀值（对这个序列来说）。</p>
</div>
</div>
<div class="section" id="id11">
<h2>14.5  模体对象自身相关的位置特异性得分矩阵<a class="headerlink" href="#id11" title="Permalink to this headline">¶</a></h2>
<p>为了更好的利用PSSMs来查找潜在的TFBSs，每个模体都同位置权重矩阵和位置特异性得分矩阵相关联。用Arnt模体来举个例子：</p>
<p>在这出现的负无穷大是由于在频率矩阵中相关项的值为0,并且我们默认使用为作为伪计数：</p>
<p>如果你更改了 <tt class="docutils literal"><span class="pre">.pseudocouts</span></tt> 属性，那位置频率矩阵和位置特异性得分矩阵都会自动地重新计算：</p>
<p>如果你想对4个核苷酸使用不同的伪计数，可以使用字典来设定4个核苷酸的 <tt class="docutils literal"><span class="pre">pseudocounts</span></tt> 。把 <tt class="docutils literal"><span class="pre">motif.pseudocounts</span></tt> 设为 <tt class="docutils literal"><span class="pre">None</span></tt> 会让伪计数重置为0的默认值。</p>
<p>位置特异性得分矩阵依赖于一个默认均一的背景分布：</p>
<p>同样，如果你更改了背景分布，位置特异性得分矩阵也会重新计算：</p>
<p>把 <tt class="docutils literal"><span class="pre">motif.backgroud</span></tt> 设为 <tt class="docutils literal"><span class="pre">None</span></tt> 将重置为均一的分布。</p>
<p>如果你把 <tt class="docutils literal"><span class="pre">motif.background</span></tt> 设为一个单一值，这个值将会被看成是GC含量：</p>
<p>应当注意到你能够在基于它计算的在背景下计算PSSM的平均值：</p>
<p>它的标准方差也是一样：</p>
<p>和它的分布：</p>
<p>请注意，每当你调用 <tt class="docutils literal"><span class="pre">motif.pwm</span></tt> 或 <tt class="docutils literal"><span class="pre">motif.pssm</span></tt> 时，位置仅重矩阵和位置特异性得分矩阵都会重新计算。如果看重速度并且你要重复用到PWM或PSSM时，你可以把他们保存成变量，像下面那样：</p>
</div>
<div class="section" id="id12">
<h2>14.6  模体比较<a class="headerlink" href="#id12" title="Permalink to this headline">¶</a></h2>
<p>当我们有多个模体的时候，我们就会想比较他们。</p>
<p>在我们开始比较之前，应当要指出模体的边界通常相当模糊。这也就是说我们需要比较不同长度的模体，因此这些比较也要涉及到相关的比对。所以我们需要考虑两个东西：</p>
<ul class="simple">
<li>模体的比对</li>
<li>比较比对后模体的相关函数</li>
</ul>
<p>为了比对模体，我们使用PSSMs的不含间隔的比对，并且用0来代替矩阵开始和结束位置缺失的列。这说明我们能够有效地利用利用背景分布来代替PSSM中缺失的列。距离函数然后可以返回模体间最小的距离，以及比对中相应的偏移量。</p>
<p>让我们举个例子，先载入另外的模体，这个模体和我们的测试模体 <tt class="docutils literal"><span class="pre">m</span></tt> 相似：</p>
<p>为了让模体能够进行相互比较，我们选择和模体 <tt class="docutils literal"><span class="pre">m</span></tt> 相同伪计数和背景值：</p>
<p>我们将用Pearson相关来比较这些模体。由于我们想要让它偏向于一个距离长度，我们实际上取1−<em>r</em> ，其中 <em>r</em> 是Pearson相关系数（PCC）：</p>
<p>这意味着模体 <tt class="docutils literal"><span class="pre">m</span></tt> 和 <tt class="docutils literal"><span class="pre">m_reb1</span></tt> 间最佳PCC可以从下面的比对中获得：</p>
<p>其中 <tt class="docutils literal"><span class="pre">b</span></tt> 代表背景分布。PCC值大概为1−0.239=0.761。</p>
</div>
<div class="section" id="de-novo">
<h2>14.7  查找 <em>De novo</em> 模体<a class="headerlink" href="#de-novo" title="Permalink to this headline">¶</a></h2>
<p>如今，Biopython对 <em>De novo</em> 模体查找只有有限的支持。也就是说，我们支持AlignAce和MEME的运行和读取。由于模体查找工具如雨后春笋般出现，所以很欢迎新的分析程序加入进来。</p>
<div class="section" id="id13">
<h3>14.7.1  MEME<a class="headerlink" href="#id13" title="Permalink to this headline">¶</a></h3>
<p>假设你用MEME以你喜欢的参数设置来跑序列，并把结果保存在文件 <tt class="docutils literal"><span class="pre">meme.out</span></tt> 。你可以用以下的命令来得到MEME输出的模体：</p>
<p>除了最想要的一系列模体外，结果中还包含了很多有用的信息，可以通过那些一目了然的属性名获得：</p>
<ul class="simple">
<li><tt class="docutils literal"><span class="pre">.alphabet</span></tt></li>
<li><tt class="docutils literal"><span class="pre">.datafile</span></tt></li>
<li><tt class="docutils literal"><span class="pre">.sequence_names</span></tt></li>
<li><tt class="docutils literal"><span class="pre">.version</span></tt></li>
<li><tt class="docutils literal"><span class="pre">.command</span></tt></li>
</ul>
<p>由MEME解析得到的模体可以像平常的模体对象（有实例）一样处理，它们也提供了一些额外的功能，可以为实例增加额外的信息。</p>
</div>
<div class="section" id="alignace">
<h3>14.7.2  AlignAce<a class="headerlink" href="#alignace" title="Permalink to this headline">¶</a></h3>
<p>我们可以用AlignACE程序实现类似的效果。假如，你把结果保存在 <tt class="docutils literal"><span class="pre">alignace.out</span></tt> 文件中。你可以用下面的代码读取结果：</p>
<p>同样，你的模体也和正常的模体对象有相同的属性：</p>
<p>事实上，你甚至可以观察到，AlignAce找到了一个和MEME非常相似的模体。下面只是MEME模体互补链的一个较长版本：</p>
<p>如果在你的机器上安装了AlignAce，你可以直接从Biopython中运行AlignAce。下面就是一个如何运行AlignAce的简单例子（其他参数可以用关键字参数来调用）：</p>
<p>由于AlignAce把所有的结果输出到标准输出，所以你可以通过读取结果的第一部分来获得模体：</p>
</div>
</div>
<div class="section" id="id14">
<h2>14.8  相关链接<a class="headerlink" href="#id14" title="Permalink to this headline">¶</a></h2>
<ul class="simple">
<li><a class="reference external" href="http://en.wikipedia.org/wiki/Sequence_motif">Sequence motif</a> in
wikipedia</li>
<li><a class="reference external" href="http://en.wikipedia.org/wiki/Position_weight_matrix">PWM</a> in
wikipedia</li>
<li><a class="reference external" href="http://en.wikipedia.org/wiki/Consensus_sequence">Consensus
sequence</a> in
wikipedia</li>
<li><a class="reference external" href="http://bio.cs.washington.edu/assessment/">Comparison of different motif finding
programs</a></li>
</ul>
</div>
<div class="section" id="bio-motif">
<h2>14.9  旧版Bio.Motif模块<a class="headerlink" href="#bio-motif" title="Permalink to this headline">¶</a></h2>
<p>这章剩下的部分将介绍Biopython 1.61版本前的 <tt class="docutils literal"><span class="pre">Bio.Motifs</span></tt> 模块，该模块取代了Biopython 1.50版本中基于两个早期Biopython模块—— <tt class="docutils literal"><span class="pre">Bio.AlignAce</span></tt> 和 <tt class="docutils literal"><span class="pre">Bio.MEME</span></tt> 的 <tt class="docutils literal"><span class="pre">Bio.Motif</span></tt> 模块。</p>
<p>为了平滑的过渡，早期的 <tt class="docutils literal"><span class="pre">Bio.Motif</span></tt> 模块将会和它的取代者 <tt class="docutils literal"><span class="pre">Bio.Motifs</span></tt> 一同维护到至少发行两个版本，并且最少一年时间。</p>
<div class="section" id="id15">
<h3>14.9.1  模体对象<a class="headerlink" href="#id15" title="Permalink to this headline">¶</a></h3>
<p>由于我们对模体分析感兴趣，不过让我们首先看看 <tt class="docutils literal"><span class="pre">Motif</span></tt> 对象。第一步要先导入模体库：</p>
<p>然后我可以开始创建我们的第一个模体对象。让我们创建一个DNA模体：</p>
<p>现在这里面什么也没有，让我往新建的模体加入一些序列：</p>
<p>现在我们有了一个完整的 <tt class="docutils literal"><span class="pre">Motif</span></tt> 实例，我们可以试着从中获取一些基本信息。先看看长度和一致序列：</p>
<p>对于DNA模体，我们还可以获得一个模体的反向互补序列：</p>
<p>我们也可以简单的调取模体的信息容量：</p>
<p>这给我们提供了模体中信息容量的比特数，知道和背景有多少不同。</p>
<p>展示模体最常用的就是PWM（位置仅重矩阵）。它概括了在模体上任意位置出现一个符号（这里指核苷酸）的概率。这个可以用 <tt class="docutils literal"><span class="pre">.pwm()</span></tt> 方法来计算：</p>
<p>模体的PWM中的概率是基于实例中的计数，但我们发现，虽然模体中没有出现G和C，可是它们的概率仍然是非0的。这主要是因为有伪计数的存在，简单地说，就是一种常用的方式来承认我们认知的不完备以及为了避免使用0进行对数运算而出现的技术问题。</p>
<p>我可以调整伪计数添加到模体对象两个属性的方式。 <tt class="docutils literal"><span class="pre">.background</span></tt> 是我们假设代表背景分布的所有字符的概率分布，是非模体序列（通常基于各自基因组的GC含量）。在模体创建的时候，就默认的设置为一个统一的分布：</p>
<p>另一个就是 <tt class="docutils literal"><span class="pre">.beta</span></tt> ，这个参数可以说明我们应该给伪计数设定为何值。默认设定为1.0。</p>
<p>所以输入伪计数的总量等于一个实例的输入总量。</p>
<p>使用背景分布和附加了伪计数的pwm，可以很容易的计算log-odd比率，告诉我们在背景下，一个来自模体特定碱基的log-odd值。我们可以使用 <tt class="docutils literal"><span class="pre">.log_odds()</span></tt> 方法：</p>
<p>在这，我们可以看到如果模体中的碱基比背景中出现频率更高，其值为正值，反之则为负值。0.0说明在背景和模体中出现的概率是相同的（如第二个位置的“T”）。</p>
<div class="section" id="id16">
<h4>14.9.1.1  模体读写<a class="headerlink" href="#id16" title="Permalink to this headline">¶</a></h4>
<p>手动从实例创建一个模体确实没什么技术含量，所以很有必要有一些读写功能来读取和写入模体。对于如何存储模体还没有一个固定的标准，但是有一些格式比其他格式更流行。这些格式的主要区别在于模体的创建是基于实例还是一些PWM矩阵。其中一个最流行的模体数据库就是 <a class="reference external" href="http://jaspar.genereg.net">JASPAR</a> ，该数据库保存了上面提到的两种类型的格式，所以让我看看我们是如何从实例中导入JASPAR模体：</p>
<p>从一个计数矩阵中导入：</p>
<p><tt class="docutils literal"><span class="pre">arnt</span></tt> 和 <tt class="docutils literal"><span class="pre">srf</span></tt> 模体可以为我们做相同的事情，但是它们使用不同的内部表现形式来展现模体。我们可以用 <tt class="docutils literal"><span class="pre">has_counts</span></tt> 和 <tt class="docutils literal"><span class="pre">has_instances</span></tt> 属性来区分它们：</p>
<p>对于模体的不同表现形式，可以用转换功能来实现相互转换：</p>
<p>在这里要注意的是 <tt class="docutils literal"><span class="pre">make_instances_from_counts()</span></tt> 方法创建的是假实例，因为按照相同的pwm能够得到许多不同的实例，所以不能够反过来重建原来的矩阵。不过这对我们利用PWM来展现模体没有什么影响，但从基于计数的模体中导出实例时要小心。</p>
<p>说到导出，让我们看看导出函数。我们可以按fasta的格式导出：</p>
<p>或者是按TRANSFAC样的矩阵格式导出（能被一些处理软件识别）</p>
<p>最后，如果我们能够联网，我们可以创建一个 <a class="reference external" href="http://weblogo.berkeley.edu">weblogo</a> ：</p>
<p>我们可以把得到的标识图以png的格式保存到特定的文件中。</p>
</div>
</div>
<div class="section" id="id17">
<h3>14.9.2  查找实例<a class="headerlink" href="#id17" title="Permalink to this headline">¶</a></h3>
<p>模体中最常用的就是在一些序列中查找实例。为了说明这章，我们将手动创建一个下面一个序列：</p>
<p>查找实例最简单的方法就是在模体中查找具体匹配的实例：</p>
<p>对于互补序列，也可以用相同的方法（为了找到互补链上的实例）：</p>
<p>It’s just as easy to look for positions, giving rise to high log-odds
scores against our motif:</p>
<p>你可能注意到阀值参数，在这里随意地设为5.0。按 <em>log</em><sub>2</sub> 来算，我们应当查找那些在模体中出现概率为背景中出现概率32倍的序列。默认的阀值是0.0,在些阀值下，会把所有比背景中出现概率大的模体实例都找出来。</p>
<p>如果你不想那么随意的选择一个阀值，你可以研究一下 <tt class="docutils literal"><span class="pre">Motif.score_distribution</span></tt> 类，它为模体提供一个相应的得分分布。由于得分的空间分布随着模体长度而成倍增长，我们正用一个近似于一个给定精度值来使计算成本变得可管理：</p>
<p>上面那个sd对象可以用来决定许多不同的阀值。</p>
<p>我们可以设定一个需要的假阳性率（找到一个由这个序列在这个背景下产生的一个模体实例的概率）：</p>
<p>或者假阴性率（找不到模体产生的实例的概率）：</p>
<p>或者一个阀值（近似），满足假阳性率和假阴性率之间的关系（fnr/fpr≃ <em>t</em>)：</p>
<p>或者一个阀值能够大体满足假阳性率和信息含量的 −<em>log</em> 值之间的相等关系（像Hertz和Stormo的Patser软件所用的一样）：</p>
<p>在我们这个例子中，当以1000比率的平衡阀值查找实例时，你可以得到一个让你获得相同结果（对于这个序列来说）的阀值：</p>
</div>
<div class="section" id="id18">
<h3>14.9.3  模体比较<a class="headerlink" href="#id18" title="Permalink to this headline">¶</a></h3>
<p>当我们有多个模体的时候，我们就会想比较他们。对此， <tt class="docutils literal"><span class="pre">Bio.Motif</span></tt> 有三种不同的方法来进行模体比较。</p>
<p>在我们开始比较之前，应当指出模体的边界通常是相当模糊的。也就是说我们经常需要比较不同长度的模体，因此这些比较涉及到相关的比对。所以我们需要考虑两个要点：</p>
<ul class="simple">
<li>模体的比对</li>
<li>比较比对后模体的相关函数</li>
</ul>
<p>在 <tt class="docutils literal"><span class="pre">Bio.Motif</span></tt> 中有三种比较方法，这些方法都是基于来源于模体比对的想法，而采用不同方式。简单来说，我们使用不含间隔的PSSMs比对，并且用0来代替矩阵同背景相比，在开始和结束位置出现缺失的列。这三种比较方法都可以解释成距离估量，但是只有一个（ <tt class="docutils literal"><span class="pre">dist——dpq</span></tt> ）满足三角不等式。这些方法都返回距离的最小值和模体相应的偏移量。</p>
<p>为了展示这些比较功能是如何实现的，让我加载和测试模体 <tt class="docutils literal"><span class="pre">m</span></tt> 相似的其他模体：</p>
<p>我们第一个展示的功能是基于Pearson相关的。因为我们想让它类似于一个距离估量，所以我们实际上取 1−<em>r</em> ，其中的 <em>r</em> 是Pearson相关系数（PCC）：</p>
<p>这意味着模体 <tt class="docutils literal"><span class="pre">m</span></tt> 和 <tt class="docutils literal"><span class="pre">Ubx</span></tt> 间最佳的PCC可以从下面的比对中获得：</p>
<p>其中 <tt class="docutils literal"><span class="pre">b</span></tt> 代表背景分布。PCC值大概为 1-0.42=0.58.如果我们尝试计算Ubx模体的互补序列：</p>
<p>我们可以发现更好的PCC值（大概为0.75），并且比对也是不同的：</p>
<p>还有两个其他的功能函数： <tt class="docutils literal"><span class="pre">dist_dpq</span></tt> ,这是基于Kullback-Leibler散度的真正度量（满足三角不等式）。</p>
<p>还有 <tt class="docutils literal"><span class="pre">dist_product</span></tt> 方法，它是基于概率的产物，这概率可以看成是两个模体独立产生两个相同实例的概率。</p>
</div>
<div class="section" id="id19">
<h3>14.9.4  <em>De novo</em> 模体查找<a class="headerlink" href="#id19" title="Permalink to this headline">¶</a></h3>
<p>目前，Biopython对 <em>de novo</em> 模体查找只有一些有限的支持。也就是说，我们只支持AlignAce和MEME的运行和读取。由于现模体查找工具发展如雨后春笋般，我们很欢迎有新的贡献者加入。</p>
<div class="section" id="id20">
<h4>14.9.4.1  MEME<a class="headerlink" href="#id20" title="Permalink to this headline">¶</a></h4>
<p>假如你以中意的参数用MEME来跑你自己的序列，并把得到的结果保存在 <tt class="docutils literal"><span class="pre">meme.out</span></tt> 文件中。你可以用以下代码读取MEME产生的文件获得那些模体：</p>
<p>除了那一系列想要的模体外，结果对象中还有很多有用的信息，可以那些一目了然的属性名来获取：</p>
<ul class="simple">
<li><tt class="docutils literal"><span class="pre">.alphabet</span></tt></li>
<li><tt class="docutils literal"><span class="pre">.datafile</span></tt></li>
<li><tt class="docutils literal"><span class="pre">.sequence_names</span></tt></li>
<li><tt class="docutils literal"><span class="pre">.version</span></tt></li>
<li><tt class="docutils literal"><span class="pre">.command</span></tt></li>
</ul>
<p>MEME解析器得到的模体可以像通常模体（含有实例）一样进行处理，它们也可以通过对实例添加附加信息而提供一些额外的功能。</p>
</div>
<div class="section" id="id21">
<h4>14.9.4.2  AlignAce<a class="headerlink" href="#id21" title="Permalink to this headline">¶</a></h4>
<p>对于AlignACE程序也可以做相同的事情。假如你把结果存储于文件 <tt class="docutils literal"><span class="pre">alignace.out</span></tt> 文件中。你可以用以下代码读取结果：</p>
<p>同样，得到的模体也和平常的模体一样：</p>
<p>事实上，你甚至可以发现AlignAce和MEME得到的模体十分相似，只不过AlignAce模体是MEME模体反向互补序列的加长版本而已：</p>
<p>如果在你的机器上安装了AlignAce，你也可以直接从Biopython中启动。下面就是一个如何启动的小例子（其他参数可以用关键字参数指定）：</p>
<p>由于AlignAce把结果打印到标准输出，因此你可以通过读取结果的第一部分来获得你想要的模体：</p>
</div>
</div>
</div>
</div>


          </div>
        </div>
      </div>
      <div class="sphinxsidebar">
        <div class="sphinxsidebarwrapper">
  <h3><a href="index.html">Table Of Contents</a></h3>
  <ul>
<li><a class="reference internal" href="#">第14章   使用Bio.motifs进行模体序列分析</a><ul>
<li><a class="reference internal" href="#id1">14.1  模体对象</a><ul>
<li><a class="reference internal" href="#id2">14.1.1  从实例中创建一个模体</a></li>
<li><a class="reference internal" href="#id3">14.1.2  读取模体</a><ul>
<li><a class="reference internal" href="#jaspar">JASPAR</a></li>
<li><a class="reference internal" href="#meme">MEME</a></li>
<li><a class="reference internal" href="#mast">MAST</a></li>
<li><a class="reference internal" href="#transfac">TRANSFAC</a></li>
</ul>
</li>
<li><a class="reference internal" href="#id4">14.1.3  模体写入</a></li>
<li><a class="reference internal" href="#id5">14.1.4  绘制序列标识图</a></li>
</ul>
</li>
<li><a class="reference internal" href="#id6">14.2  位置权重矩阵</a></li>
<li><a class="reference internal" href="#id7">14.3  位置特异性得分矩阵</a></li>
<li><a class="reference internal" href="#id8">14.4  搜索实例</a><ul>
<li><a class="reference internal" href="#id9">14.4.1  搜索准确匹配实例</a></li>
<li><a class="reference internal" href="#pssm">14.4.2  用PSSM得分搜索匹配实例</a></li>
<li><a class="reference internal" href="#id10">14.4.3  选择一个得分阀值</a></li>
</ul>
</li>
<li><a class="reference internal" href="#id11">14.5  模体对象自身相关的位置特异性得分矩阵</a></li>
<li><a class="reference internal" href="#id12">14.6  模体比较</a></li>
<li><a class="reference internal" href="#de-novo">14.7  查找 <em>De novo</em> 模体</a><ul>
<li><a class="reference internal" href="#id13">14.7.1  MEME</a></li>
<li><a class="reference internal" href="#alignace">14.7.2  AlignAce</a></li>
</ul>
</li>
<li><a class="reference internal" href="#id14">14.8  相关链接</a></li>
<li><a class="reference internal" href="#bio-motif">14.9  旧版Bio.Motif模块</a><ul>
<li><a class="reference internal" href="#id15">14.9.1  模体对象</a><ul>
<li><a class="reference internal" href="#id16">14.9.1.1  模体读写</a></li>
</ul>
</li>
<li><a class="reference internal" href="#id17">14.9.2  查找实例</a></li>
<li><a class="reference internal" href="#id18">14.9.3  模体比较</a></li>
<li><a class="reference internal" href="#id19">14.9.4  <em>De novo</em> 模体查找</a><ul>
<li><a class="reference internal" href="#id20">14.9.4.1  MEME</a></li>
<li><a class="reference internal" href="#id21">14.9.4.2  AlignAce</a></li>
</ul>
</li>
</ul>
</li>
</ul>
</li>
</ul>

  <h4>Previous topic</h4>
  <p class="topless"><a href="chr13.html"
                        title="previous chapter">Chapter 13  Phylogenetics with Bio.Phylo</a></p>
  <h4>Next topic</h4>
  <p class="topless"><a href="chr15.html"
                        title="next chapter">Chapter 15  Cluster analysis</a></p>
  <h3>This Page</h3>
  <ul class="this-page-menu">
    <li><a href="_sources/chr14.txt"
           rel="nofollow">Show Source</a></li>
  </ul>
<div id="searchbox" style="display: none">
  <h3>Quick search</h3>
    <form class="search" action="search.html" method="get">
      <input type="text" name="q" />
      <input type="submit" value="Go" />
      <input type="hidden" name="check_keywords" value="yes" />
      <input type="hidden" name="area" value="default" />
    </form>
    <p class="searchtip" style="font-size: 90%">
    Enter search terms or a module, class or function name.
    </p>
</div>
<script type="text/javascript">$('#searchbox').show(0);</script>
        </div>
      </div>
      <div class="clearer"></div>
    </div>
    <div class="related">
      <h3>Navigation</h3>
      <ul>
        <li class="right" style="margin-right: 10px">
          <a href="genindex.html" title="General Index"
             >index</a></li>
        <li class="right" >
          <a href="chr15.html" title="Chapter 15 Cluster analysis"
             >next</a> |</li>
        <li class="right" >
          <a href="chr13.html" title="Chapter 13 Phylogenetics with Bio.Phylo"
             >previous</a> |</li>
        <li><a href="index.html">Biopython_en 1.0 documentation</a> &raquo;</li> 
      </ul>
    </div>
    <div class="footer">
        &copy; Copyright 2013, Biopython.
      Created using <a href="http://sphinx-doc.org/">Sphinx</a> 1.2b1.
    </div>
  </body>
</html>